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ABSTRACT 


^""^A  new  technique  for  numerically  solving  the  reduced 
wave  equation  on  exterior  domains  is  presented.  The  method 
is  basically  a  relaxation  scheme  which  exploits  the  limiting 
amplitude  principle.  A  modified  boundary  condition  at  "infinity'.' 
is  also  given.  The  technique  is  tested  on  several  model 
problems:  the  scattering  of  a  plane  wave  off  a  metal  cylinder, 

a  metal  strip,  a  Helmholtz  resonator,  an  inhomogeneous 
cylinder  (lens  ),  and  a  nonlinear  plasma  column., 

The  results  are  in  good  qualitative  agreement  with  previously 
calculated  values.  In  particular,  the  numerical  solutions 
exhibit  the  correct  refractive  and  diffractive  effects  at 
moderate  frequencies.-,. — 

1.  Introduction 

It  is  well  known  that  for  dissipative  linear  ordinary 
differential  equations  with  a  forcing  term  of  period  \  that  the 
transients  die  out  and  the  solution  tends  as  t  -*■  00  to  solutions 
of  period  A.  The  same  is  true  of  many  hyperbolic  equations.  In 
particular  solutions  of  the  wave  equation  in  infinite  domains 
with  appropriate  boundary  conditions  and  with  the  forcing  term 
felu)t  tend,  for  large  times,  to  solutions  of  the  Helmholtz 
equation  Au  +  co  u  =  f.  This  is  known  as  the  limiting  amplitude 
principle  [10], 
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The  purpose  of  this  paper  is  to  show  that  a  varied  form 
of  a  limiting  amplitude  principle  can  be  used  to  solve 
numerically,  in  a  short  time  and  for  a  rather  wide  variety 
of  geometries,  the  Helmholtz  equation  in  the  exterior  of  an 
obstacle  with  a  variable  index  of  refraction.  Furthermore 
this  can  be  done  at  such  high  frequencies  that  a  great 
deal  of  geometric  optical  behavior  can  be  confirmed  even 
on  a  relatively  coarse  mesh. 

The  Helmholtz  equation  with  variable  index  of  refraction 
is 

2 

Au  +  a)  nu  =  0 


where  the  index  of  refraction  n  =  n(x)  is  a  function  of  the 
space  variables  (x).  We  also  consider  situations  where 
n  =  n(x,|u|)  which  corresponds  to  certain  models  of  laser 
beam  propagation.  At  large  distances  we  assume  n  •*  constan 
which  may  be  scaled  to  1.  Then  the  parameter  to  has  the 
dimension  of  [Length]  x  and  the  appropriate  dimensionless  r ■<* 

constant  is  ua  where  a  is  a  characteristic  length  of  the 


scatterer.  It  may  range  from  zero  to  infinity. 


The  field  u  is  given  as  a  plane  wave  e 
scattered  field  ug; 


1U)X 


and  its 

/ 


(l.D 


110X  , 

u  =  e  +  u. 


ft 


However,  the  source  could  equally  well  be  located  at  a  finite 


*wm 


point  in  the  plane.  The  scattered  field  satisfies  the  outward 
radiation  condition 


(1.2) 


8u  u 

s  ,  s 

TT -  ~  ltOU  +  7T- 

3r  s  2r 


as  x  -*■ 


This  problem  may  be  studied  using  geometrical  optics 
for  high  frequency  and  expansion  in  w  for  low  frequency. 

Our  original  objective  was  to  deal  with  the  range  of  frequency 
where  neither  of  these  approximations  is  good.  We  found,  in 
addition,  however,  that  many  of  the  features  of  geometrical 
optics  emerge  from  the  computations  at  moderate  frequencies. 
Primarily  though,  we  are  concerned  with  perturbations  in  the 
index  of  refraction  or  with  disturbing  objects  where  the 
characteristic  length  is  not  large  compared  to  the  wave  length. 

A  wave  at  frequency  w  has  the  wave  length  2nVa)  and  theo¬ 
retically  for  geometrical  optics  one  needs  2tt/o)  <<  a.  Our  ob¬ 
servations  confirm  that  in  many  applications  wa  =  5  or  10 
displays  the  significant  features  of  geometrical  optics  and 
diffraction  theory. 

In  this  paper  we  restrict  ourselves  to  two  dimensional 
examples.  In  a  second  paper  we  plan  to  demonstrate  results 
in  axi-symmetric  and  possibly  three  dimensional  geometries. 

A  particular  advantage  in  2D  is  that  every  simply-connected 
object  can  be  studied  by  noting  that  it  can  be  mapped  con¬ 
formally  onto  the  exterior  of  a  circle  and  the  transformation 
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induces  a  .new  Helmholtz  equation  for  the  solution  with  a  new 
index  of  refraction.  In  fact,  however,  we  have  also  computed 
objects  that  are  not  simply-connected  but  have  rectangular 
boundaries  in  polar  coordinates. 

The  principal  results  of  interest  are: 

(i)  Confirmation  of  the  scattering  cross  section  for  a 

cylinder  of  radius  one  with  a  Dirichlet  boundary  con¬ 
dition.  Results  are  compared  with  those  presented  by 
Bowman,  Senior,  and  Uslenghi  [  i  D .  The  error  is  less 
than  1(10)  %  at  the  frequency  to  =  5  0.0)  .  The 

computations  also  show  that  the  position  of  the  shadow 
edge,  as  given  by  geometrical  optics,  is  within  one 
or  two  mesh  points. 

(ii)  The  diffraction  by  an  infinite  strip  of  half-width 

2  at  various  angles  to  an  incident  plane  wave.  There 
is  good  correlation  with  [  1  ]  particularly  at  low 
frequencies  and  at  an  angle  of  45°  and  with  geometrical 
optics  [  1  ]  . 

(iii)  The  computations  of  the  field  including  the  cross- 
section  of  a  Helmholtz  resonator  of  radius  1  with 
various  apertures  and  for  plane  waves  at  various  angles 
of  incidence. 


(iv)  The  refraction  of  a  plane  wave  by  a  lens  of  either 
constant  or  variable  index  of  refraction. 


In  particular  we  include  lenses  that  produce  focusing  and 
caustics  such  as  the  Luneberg  lens. 

(v)  The  refraction  by  a  model  of  an  overdense  plasma. 

The  method  of  computation  is  a  modification  of  the 
relaxation  method  used  in  [  2  ]  and  is  described  in  Section  2. 
To  reduce  the  number  of  mesh  points  involved  a  modified 
radiation  condition  corresponding  to  (2)  is  imposed  at  a 
finite  radius  (see  Section  3).  The  difference  scheme  and 
how  to  deal  with  the  artificial  singularity  at  the  origin 
created  by  using  polar  coordinates  is  described  in  Section  4. 
In  Section  5,  we  describe  in  detail  the  particular  examples 
computed.  In  each  case  we  give  the  running  times  on  the  CDC 
6600  required  to  achieve  convergence.  The  effects  of  apply¬ 
ing  the  modified  radiation  condition  at  different  distances 
are  also  discussed. 

In  Section  6  we  discuss  the  limitations  of  this  kind 
of  iterative  method  and  other  possible  applications. 

The  authors  are  grateful  to  F.  Tappert  for  much  advice 
and  many  suggestions.  The  work  of  A.  Bayliss  and  E.  Turkel 
[  3  ]  on  similar  problems  with  obstacles  treats  the  radiation 
condition  using  the  same  ideas.  Turkel  pointed  out  to  the 
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authors  that  the  modified  radiation  condition  is  well  posed 
in  the  sense  of  Kreiss-Lobatchinskii  [  11  3 . 

2.  The  Iteration  Method. 

We  consider  a  time  dependent  equation 

(2.1)  (QU)t  =  AU  +  mo2U 

where  Q  is  a  first  order  operator  in  the  space  variables, 

Q  =  2a- 7  +  b.  The  coefficient  vector  a  and  the  scalar  b  are 
chosen  so  that  the  solutions  of  (2.1)  of  the  form  .plane  wave 
or  source  plus  outgoing  scattered  wave  will  approach  a  time 
independent  state.  This  will  be  the  desired  solution. 

This  method  was  used  in  [23.  However,  (2.1)  was  used 
directly,  i.e.  in  characteristic  form.  Data  was  given  on  a 
characteristic  surface.  The  convergence  was  therefore  very 
slow  because  a  very  small  time  step  was  required  for 
stability  with  the  space  differences  used. 

In  the  present  method  we  transform  (2.1)  to  a  Cauchy 
problem  by  a  change  of  variables  given  by 

dt '  =  dt  +  a-dx 

(2.2) 

dx'  =  dx. 

We  obtain  in  the  new  variables 

(2.3)  1  a | 2^tt  =  AU  +  Ut(V*a-b)  +  nu2U 


where  we  have  dropped  the  primes.  We  now  solve  an  initial 
value  problem  which  by  our  choice  of  a  and  b  will  converge 
to  the  steady  state.  For  example,  the  convergence  has  been 
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improved  by  more  than  an  order  of  magnitude  in  the  case 


(2.4) 


a  =  x/|x|,  b  =-2iw  +  V-a. 


In  two  dimensions  with  b  given  by  (2.4)  we  introduce 
in  (2.3)  the  change  of  variable: 


(2.5) 


U  =  eiwx  +  ^ 


with  |x|  =  r.  Equation  (2.3)  reduces  in  polar  coordinates  to 

(2.6)  W.  .  =  W  +r~2 tWQa+%W]  t  a)2 (n-1  )W  +  w2/r(n-l)elco(x-t) 
tt  rr  v u 


which  for  n  =  1  is  the  wave  equation.  However,  in  general 
it  has  fixed  characteristics  that  are  independent  of  n.  Thus 
if  n  =  1  we  are  using  the  standard  limiting  amplitude  principle 
but  since  in  some  of  our  examples  n  changes  dramatically  it 
is  a  great  advantage  that  the  characteristics  are  fixed.  For 
n  i  1  we  are  solving  a  wave  equation  with  potential. 

From  [  2  ]  we  have  the  growth  restriction  to  obtain  decay 


(2.7) 


■gp[r(n-l)]  >  0 


If  this  condition  is  not  satisfied  we  are  dealing  with  a  po¬ 
tential  which  has  bound  states  that  give  rise  to  exponen¬ 
tially  growing  solutions  of  (2.6).  This  showed  up  numerically 
in  a  very  dramatic  fashion  when  we  tried  lenses  with  n  >  1. 

The  inequality  (2.7)  was  proved  by  J.  Weidmann  [4]  as  a  necessary 
condition  to  prevent  these  growing  modes. 
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In  order  to  handle  this  situation  we  note  that  instead 
of  (2.6)  we  could  use 

(2.8)  ni^tt  =  ^rr  +  r  +  (n^-Du^W  +  /r (n-1  )e^w ^x-t ^ 

where  n^  +  ^  =  n  +  1.  If  the  limiting  amplitude  principle 
holds,  i.e.  the  potential  (n2“l)w  doe:;  ~>ot  give  rise  to 
growing  modes  and  n^  >  0 ,  then  the  solutions  of  (2.8)  also 
approach  a  time  harmonic  steady  state  and  the  space  dependent 
coefficient  satisfies  the  desired  Helmholtz  equation. 

If  we  make  n^  >  1 

in  such  a  way  that  [r(n2~l)]r  s  0,  then  (2.8)  can  be  used  for 
the  iteration.  We  have  applied  this  notion  only  for  n^  =  n  >  1, 
n2  =  1  and  n  <  4  inside  a  circle.  Note  that  the  Courant- 
Friedrichs-Lewy  condition  for  the  time  step  continues  to  be 
satisfied  over  n^  >  1  if  it  holds  for  n^  =  1 .  On  the  other 
hand  if  there  are  trapped  rays  produced  by  the  artificial 
index  n-^  then  the  approach  to  steady  state  may  be  exceedingly 
slow  as  the  higher  frequencies  have  poorer  exponential  decay. 

3.  Improving  the  Radiation  Condition  and  Finding  the 
Scattering  Cross  Section. 

Returning  to  the  original  radiation  condition  (1.2)  we 
note  that  if  a  solution  satisfies  the  radiation  condition, 
then  by  using  the  fundamental  solution  representation  for 
the  reduced  wave  equation,  we  have  with  us  =  Wr  ^ ^ ^ elu>r 

the  expansion 
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(3.1) 


W  =  E  A.r 

5=o  ^ 


in  any  number  of  dimensions.  Here  A^  is  a  function  of  w  and 
the  angular  variables.  The  differential  equation  for  W  at 
large  distances  is 


(3.2) 


2iWr  =  -Wrr  -  ±5-  LW 
r 


where  L  is  a  differential  operator  acting  on  the  angular 
variables  (the  Laplace-Beltrami  operator).  We  have  rescaled 
r  so  that  w  =  1.  On  substituting  the  expansion  for  W  we 
find  the  recursion  formula 


(3.3)  An+]_  =  P(n+l,L)An,  n  >  0 

where  P(n+1,L)  =  -  yi  (n+1)  *  (L+n(n+l)}.  Note  that  PQ,L)  =  ~ 
Next  we  set 


n 

P(n,L)  =  n  P(j , L) ,  P(0,L)  =  1 

5  =  1 


and  obtain 

00 

(3.4)  W  =  E  P(n,L)r"nAn 

0  U 

and 

00  , 

Wr  =  -E  P(n,L)nr"n"iA0 

which  we  invert  to  find 


A 


0 


-(E  nP(n,L)r-n"1)"1W 
1  r 
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arid  thus 

oo  oo 

(3.5)  W  =  -(£  n(n  +  l )  P(n ,  L)r  n  ^)(Z  nP(n,L)r  n”^)""^W  . 

1  1  r 

The  boundary  condition  on  W  is  then  found  by  substituting  in 
the  differential  equation  (3.2) 

00 

(3.6)  2iW  =  -  LW  +  i[2P(l,L)  +  E  n  (n+1 )  P(n ,  L)  r'  r‘+1  ]  • 

r  ^  n=2 

OO 

[P(l ,  L)  +  E  nP(n,L)r"n+1]“1v-7 
n=2  r 

and  solving  for  W  as  a  function  of  LW. 

We  simply  expand  the  formula  in  powers  of  r  but  it  may 

be  better  to  use  the  approach  of  Enquist  and  Majda  [5]  and  use 

other  representations  of  the  operators.  What  is  needed  is  a 

OO 

good  representation  of  E  nP(n,L)r  n+'*‘.  The  first  approximation 

n=2 

is  Wr  =  0  and  from  (3.6)  the  next  approximation  is 

(.3.7)  (2i  -  1)W  =  +  0(r"4). 

r 

Note  that  as  to  -»■  «,  LW  00  in  general  and  there  are  problems  of 
convergence . 

In  transforming  to  the  Cauchy  problem  the  derivative 
W^  ->  W^  +  and  the  right  hand  side  is  unchanged.  Thus  the 
radiation  condition  which  we  used  was 


(3.8) 


Wr  + 


wt  = 


LW 


2o>r  (i  -  — : 

wr 


where,  for  the  two-dimensional  problem  L 


1/4. 


38 


■f  . 
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To  obtain  the  far  field, 
section  Ag,  we  invert  (3.4) 

(3.9)  Ag  =  tl  - 

which  is  accurate  to  0(—i — . 

to  r 


that  is  the  scattering  cross 
and  obtain 

P(l,L)y.T 

tor 


4.  The  Difference  Scheme. 

To  solve  the  time-dependent  differential  equation  (1.7) 
imposing  the  far  field  condition  (3.9)  we  have  used  a  standard 
backward  difference  scheme  for  the  initial  value  problem  of  a 
second  order  hyperbolic  equation  with  second  order  accuracy. 
Let  the  solution  to  the  difference  scheme  be  W(j,m,n)  where 
(j ,m,n)  are  evaluated  on  a  grid 

(4.1)  r  =  jir  +  rQ ,  9  =  mA9,  t  =  nAt, 

with 

0  <  j  <  N,  0  s  m  <  M 


Then  for  an  interior  point , 

(4.2)  W(j,m,n+1)  =  T( j ±1 , j ,m±l ,m,n ,n-l ) 

is  determined  from  the  values  of  W  at  (j±l,m±l,n),  (j,m,n), 
(j,m,n-l).  On  the  outer  boundary  r  =  NAr  +  r^  we  use  the 
differenced  form  of  (3.9) 

(4.3)  {W(N,m,n+l  )-W(N  ,m,n-l )  }  +  yjp{W(N+l  ,m,n)-W(N-l  ,m,n 

=  S(N,m,n) 

where  S(N,m,n)  denotes  the  difierence  approximation  for 
the  right  hand  side,  centered  differences  being  used  in  L. 
The  value  of  W  at  j  =  N+l,  m,n  has  to  be  eliminated  by  using 
the  difference  equation  (4.2)  and  thus  one  obtains 

(4.4)  W(N,m,n+l)  =  B(N ,m ,m±l ,n ,n-l ) . 

It  still  remains  to  apply  the  Dirichlet  condition  if 
there  is  a  disturbing  object  or  to  deal  with  the  singularity 
at  the  origin  which  is  artificially  created  by  the  use  of 
polar  coordinates. 

(a)  To  apply  the  Dirichlet  condition  (u=0)  we  simply 
use  the  given  values  from  W  =  -/r  cxp[iu>(x-t )  ]  at  the  mesh 
point  j  =  0  for  a  circle  or  at  any  other  boundary  point.  To 
simplify  the  computation  we  have  used  only  boundaries  con¬ 
sisting  of  sections  9  =  constant  or  r  =  constant. 


(b)  To  deal  with  the  origin  we  note  two  difficulties. 
First  the  c ourant-Friedrichs-Lewy  condition  limits  the 
smallest  radius,  i.e.,  we  must  have  rA0  >  At  for  all  points 
in  the  computation.  Secondly,  there  is  a  singularity  from 
the  term  \  r  W.  However,  we  do  have  W  ~  /r .  Our  approach 
is  to  keep  doubling  the  mesh  size  A9  as  r->  0.  Thus  in  these 
problems  for  r  <  rg  we  use  the  mesh 

(4.5)  r=  j'Ar,  0  =  ms(j'),  t  =  nAt 

where  s(j')  is  chosen  so  that  ^A0  =  j'Ars(j')  >  At  and  so 
that  for  each  j'  where  the  mesh  size  s(j')  changes  the 
0-mesh  is  half  as  dense  as  before.  The  behavior  of  the 
equation  at  the  origin  itself  has  to  be  taken  into  account. 
Various  integral  forms  of  the  equation  could  be  used  but 
it  is  sufficiently  effective  to  use  the  equation  for  U  at 
r  =  0  using  for  AU  the  values  of  $  =  W//r  at  0  =  0 ,  ±tt/2,  it 
and  r  =  2Ar.  By  the  standard  difference  formula  along  with 
<J(n)  for  the  value  of  If  at  r  =  0  and  time  t  =  nAt.  This 
yields,  for  the  symmetric  case, 

(<Kn+l)-2<Mn)+<Mn-l))/(At)2 

=  (W(2,m-pn)  +  2W(2,m2,n)+W(2,mg,n)-4<J>(n))/(2Ar)^2 
+  u)2n(0  )$(n)  +  w2(n(0)-l)e 

Here  m^ ,  m2,  m^  correspond  to  0  =  0,j,tt.  Note  n(0)  means 
the  value  of  the  index  of  refraction  at  the  origin.  In  the 
non-symmetric  case  there  is  a  similar  formula. 
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Finally  we  obtain  W(l,m,n+1)  by  interpolating  the  values 

of  <5^n+1^  and  W (  2  , m , n  + 1 )  /  /2  A  r  to  obtain  W(  1  ,m,n+ 1)  / /Ar .  The 
2 

error  is  0(Ar  )  but  turns  out  to  be  larger  than  in  the  rest 
of  the  computation. 

In  closing  this  section  we  shall  describe  our  numerical 
method  for  determining  when  V/  has  reached  its  time  harmonic 
steady  state.  Recall  that  the  scattered  field,  ug  ,  is  given 
by 


us  =  [We^tlJ"t]//r. 


The  bracketted  term  must  approach  Wela,r  as  t  °°.  Thus  f  "r 
large  time  |W|  becomes  independent  of  t.  We  terminate  our 
computations  when 


(4.6) 


max 
0<i<N 
0<j  <M 


|W(n+l,i,j) |-[W(n,i,j) | 


<  e 


for  some  prescribed  e  >  0. 


5.  Results  for  Model  Problems. 


(i)  The  metal  cylinder. 

This  is  a  benchmark  problem  to  test  the  method.  Solu¬ 
tions  have  been  determined  by  separation  of  variables,  integral 
equations,  and  by  geometrical  optics  [see  reference  1  for  a 
fairly  comprehensive  bibliography].  In  [1]  the  "scaled  cross 

section"  Aoi  S(0)  is  presented  graphically  for  various  frequencies. 
2 

Here  S(6)  is  the  amplitude  of  the  outgoing  cylindrical  wave  and 
0  is  the  polar  angle  with  0  =  tt  the  direction  of  the  incident 
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plane  wave.  These  results  were  carefully  converted  to  tabular 
form  for  the  case  oia  =  5 .  They  are  shown  in  column  C  of  Table  1 
while  the  results  of  our  present  calculation  are  listed  in  column 
B.  The  areement  is  excellent.  The  relevant  parameters  used  in 
the  computation  are  wa  =  5,  Ar  =  0.1,  A9  =  u/40.  At  =  0.05,  and 
e  =  0.01. 

For  the  sake  of  comparison  the  results  of  [2]  are  listed  in 
column  A.  The  gain  over  our  old  method  is  in  the  running  time 
which  is  now  approximately  1  minute  on  the  CDC  6600.  The  old 
calculations  took  roughly  15  minutes  on  the  same  machine.  Both 
methods  used  the  same  amount  of  core,  13 8K. 

The  effect  of  applying  the  modified  radiation  condition  (3.9) 
at  various  distances  to  gain  in  mesh  refinement  was  done  by  in¬ 
creasing  the  radius  of  the  cylinder  and  applying  (3.9)  at  a  fixed 
radius.  The  relevant  parameter  is  ua  where  a  is  the  radius  of  the 
object  [see  Table  1, columns  D-F].  This  does  not,  of  course, 
refine  the  6  mesh. 

Since  our  numerical  method  gives  the  total  field  at  each 
grid  point,  the  cross  sections  given  in  Table  1  represent  a  small 
fraction  of  the  generated  information.  Instead  of  just  listing 
these  numbers  an  alternate  method  was  devised  to  visually  convey 
the  results.  First,  the  polar  output  was  converted  into  a 
rectangular  grid  of  numbers  using  straightforward  interpolation. 
(This  unfortunately  introduces  errors  which  tend  to  smear  out  the 
optical  features  that  will  be  described  shortly.)  Then  seven 
weights  of  shade  were  used  with  the  darkest  color  corresponding 


Co  the  smallest  total  field  while  the  lightest  gives  the  largest 
field.  The  amplitude  ranges  are  ( 0 , 2  )  ,  ( 2  ,  6 )  ,  ( 6 , 8 )  ,  ( 8 , 10 )  ,  ( 10 , 14 ) , 
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(l'l, 17),  and  (17,™).  This  process  gives  the  interference  pattern 
shown  in  Figure  1  when  wa  =  10.  The  light  regions  correspond 
to  constructive  interference,  the  dark  ones  to  destructive. 

Since  the  incident  wave  enters  from  \  ®  (6  =  it  or  from  the 

left  of  the  figure),  the  total  field  is  symmetric  about  the  x 
axis  and  only  half  the  pattern  is  shown  in  Figure  1.  The  dark 
semicircle  is  the  metal  cylinder.  The  wave  patterns  are  readily 
seen  in  this  picture.  Moreover,  tie  shadow  cast  by  the  cylinder 
is  quite  apparent.  The  width  of  the  transition  re  '  Ion  which 
connects  the  deep  shadow  and  the  illuminated  portions  of  the  plane 
is  exagerated  by  the  interpolation  process  mentioned  above. 


(ii)  Infinite  Strip.  (Dirichlet  Case) 

The  diffraction  by  an  infinite  strip  of  half-width 
2,  {x  =  0,  |  y  |  <  2}  was  computed  at  the  frequency  5,  u>a  =  10  at 
the  angles  of  attack  a  =  TO0  (head-on),  45°,  and  90°  (on  edge). 
Here  a  is  the  angle  between  the  incident  plane  wave  and  the 
positive  x  axis.  The  wave  again  is  approaching  from  0  =  it  . 

The  mesh  size  is  tt/20  in  9  and  0.1  in  r.  A  variable  mesh  was 
used  near  the  origin. 

Since  the  interference  pattern  in  the  x  -  y  coordinates 
is  smeared  by  interpolation  in  this  and  in  other  cases  with  sharp 
boundaries,  it  is  omitted  here.  However,  the  results  for  a  =  0 
are  presented  in  polar  form  in  Figure  2.  The  numbers  printed  at 
the  mesh  points  are  ten  times  the  total  field.  Thus  the  dark 
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regions  now  correspond  to  constructive  interference,  the  light 
to  destructive.  In  this  picture  one  sees  the  shadow  at  0,  2v , 
the  standing  wave  in  the  illuminated  region  along  Q  =  u  ,  the 
singular  corners  at  0  =  ~  and  the  shadow  boundaries  on 

y  =  t  2  =  r  sin  6. 

For  the  sake  of  brevity,  we  have  not  included  the  polar 
output  for  the  cases  a  =  45°  and  90°.  Rather  a  few  words  describing 
these  results  will  be  offered  instead.  A  turn  through  45°(=o) 
destroys  the  symmetry  and  slightly  distrots  the  waves  in  the 
illuminated  region.  The  shadow  is  shifted  45°.  When  the  wave 
attacks  the  strip  on  edge  (a  =  90°)  there  is  no  tnadow.  However, 
in  the  forward  scattered  direction  the  total  field  is  cut  in 
half  as  it  should  be;  see  [6].  The  field  is  again  symmetric. 

The  cross  sections  are  presented  in  Table  2.  They  check 
well  at  45°  and  qualitatively  at  0°  with  those  given  in  [1].  The 
deviation  between  our  results  and  those  given  by  [1]  is  probably 
due  to  the  fact  that  the  derivatives  of  the  field  became  singular 
at  the  strip's  edges.  The  coarseness  of  our  mesh  masks  this  problem 

The  running  time  and  core  requirement  for  this  problem  were 
roughly  the  same  as  those  needed  for  the  metal  cylinder  problem. 

(iii)  The  Helmholtz  Resonator 

A  cylindrical  Helmholtz  resonator  was  placed  in  a  plane  wave 
at  different  angles  of  attack.  The  resonator  is  an  infinitely 
thin  cylinder  of  radius  two  with  a  strip  aperture  centered  on  the 
negative  x  axis.  The  angle  6  which  subtends  the  aperture  was  taken 
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at  different  sizes.  When  6  =  0  the  resonator  is  a  metal  cylinder 
which  was  used  for  comparison.  (The  frequency  was  fixed  at  five.) 

The  cross  sections  arc  shown  in  Table  3  where  a  is  again  the  angle 
between  the  incident  plane  wave  and  the  positive  x  axis.  The 
other  relevant  parameters  are  At  =  0.05,  Ar  =  0.1,  and  AO  =  tt/20. 
There  is  no  data  for  comparison.  In  iterating  on  i'.-*  closed 
resonator  the  method  generated  some  eigenmodes  in  the  interior. 

These  exist  as  slowly  damped  modes  as  the  aperture  opens. 

At  an  angle  of  attack  of  a  =  0°  and  an  aperture  of  36°  the 
exterior  field  is  very  similar  to  that  generated  by  a  metal  cylinder. 
The  interior  field  contains  a  considerable  amount  of  energy.  Th^ 
amplitudes  focus  on  the  axis  in  two  places.  The  strongest  (maximum 
amplitude  5)  at  x  =  -.5  is  a  portion  of  caustic  caused  by  the  second 
reflections  off  the  interior.  The  weaker  focus  (maximum  amplitude  3) 
at  x  =  +.5  probably  corresponds  to  the  peak  of  a  nephroid-like 
caustic  formed  by  the  first  reflections.  The  difference  in  strength 
is  probably  due  to  constructive  interference. 

When  the  wave  strikes  the  resonator  from  directly  behind 
(angle  of  attack  =  180°)  it  acts  like  a  metal  cylinder.  The  only 
energy  inside  is  diffractive  and  weak. 

At  an  aperture  of  180°  and  an  angle  of  attack  of  0°  there 
should  be  a  nephroid- like  caustic  but  the  edge  diffraction  obli¬ 
terates  this  feature.  There  is  marked  focusing  as  expected, 
for  -2  s  x  s  -1.  The  shadow  is  not  so  sharp  as  with  a  metal 
cylinder.  The  aperture  edges  act  like  slits  and  smear  the 
geometrical  effects. 

The  core  requirements  for  these  problems  were  the  same 
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as  those  needed  for  the  various  metal  cylinders.  However, 
the  running  time  depends  upon  the  aperture  size  which 
determines  the  complex  eigenfrequencies  of  the  resonator. 

Several  eigenfrequencies  for  various  resonators  are  computed 
numerically  and  presented  in  [  7  ].  For  an  aperture  of  36° 
and  800  iterations  we  could  satisfy  (4.6)  only  for  e  >  .1. 
However,  when  the  aperture  was  180°  we  could  satisfy  (4.6) 
with  e  =  .05  at  400  iterations.  These  two  numerical  examples 
show  how  the  iteration  scheme  depends  upon  the  eigenfrequency 
with  the  smallest  imaginary  part.  In  the  first  case  the  imagi¬ 
nary  part  is  roughly  0.015  while  in  the  second  case  it  is  about  . 
If  we  assume  that  the  transients  decay  like  Exp(-.015t)  in 
the  first  situation,  then  t  ~  200  would  make  this  factor 
0(1/100).  Thus  for  At  =  .05  (the  number  used  in  our  program) 
we  would  need  in  the  neighborhood  of  5000  iterations  to 
obtain  convergence.  The  same  crude  argument  shows  that  about 
310  iterations  are  needed  when  the  aperture  is  180°. 

All  of  these  problems  could  be  solved  with  the  addition 
of  a  variable  index  of  refraction  depending  on  all  variables 
including  the  amplitude  of  the  total  field  provided  the 
focusing  effects  are  weak. 

(iv)  Lens  with  variable  index  of  refraction. 

For  a  lens  of  constant  index  of  refraction  less  than  1 
we  used  (2.6),  i.e.  (2.8)  with  n^  =  1,  nj  =  n.  There  are  no 
particular  difficulties  and  the  fields  are  qualitatively 


correct.  The  basic  lens  was  of  radius  a  <  3  and  the  frequency 

co  c  5  so  that  the  relevant  parameter  wa  <  15.  For  constant  n 
less  than  1  in  the  lens  there  is  rapid  convergence  and  typical 
lens  patterns  emerge.  In  Figure  3  the  interference  pattern  (in 
the  x-y  plane)  is  shown  for  a  lens  with  a  =  2,  w  =  5,  and  n  =  0.4. 

For  1  <  n  <  2  there  are  focusing  effects;  e.g.  when  n  =  2  the  maximum 
amplitude  is  2.6  and  occurs  on  the  axis  0=0.  The  full  limiting 
amplitude  principle  with  n^  =  n  and  n2  =  1  works  well  for  these 
cases.  However,  =  n  cannot  be  taken  too  large.  A  rescaling  of 
time  in  (2.S)  would  introduce  a  large  effective  frequency  into  the 
exponential  term.  This  would  generate  a  large  truncation  error  for 
a  fixed  At  and  cause  numerical  instabilities.  Numerical  experiments 
confirmed  this  observation  for  n-^  =  n  >  8  and  At  =  .05.  On  the 
other  hand,  taking  n^  =  1,  n^  -  n  brings  immediate  instability  into 
the  computation. 

More  interesting  effects  occur  if  n  =  n(r)  or  n  =  n(r,0).  For 
example,  geometrical  optics  predicts  [see  reference  8]  that  for 
n  =  r?/9,  r<3,n=l,  r>3  that  there  will  be  focusing  at  (3,^). 

At  u>  =  5  we  obtained  an  amplification  of  2.5  in  the  total  field  at 
that  point.  Also  there  is  a  geometrical  shadow  on  0  =  This 

shadow  (for  a  finite  u>  )  is  not  sharp  because  a  single  ray  on 
0  =  0 ,  0  =  n  passes  through  it. 

2 

A  second  case,  a  Luneberg  lens  with  n  =  2  -  r  /9  for  r  s  3, 
n  =  1  for  r  2  3  focused  at  (3,0)  with  a  magnitude  of  3.  This 
point  is  a  focus  for  the  geometrical  optics  rays. 

2  2  2 

The  general  case  n  =  n(r,0)  was  tried  with  n  =  (r  -  rQ)/(9  -  rQ 
for  r  s  3,  n  =  1  for  r  2  3,  where  rQ  =  2  -  0.5  cos 2 ( 0 / 2 > .  The 
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particular  effects  are  not  of  interest.  Obviously  the  oscillations 
in  n  cannot  be  too  large  without  destroying  accuracy.  What  is 
important  is  that  this  problem  cannot  be  solved  by  separation  of 
variables.  Geometrical  optics  is  complicated  and  an  integral  method 
would  involve  inverting  a  kernel  in  four  variables.  No  difficulties 
are  encountered  here.  There  is  some  increa  ,  j  in  memory  since  n 
depends  on  two  variables.  The  running  time  was  again  about  a 
minute.  Condition  (4.6)  was  satisfied  for  c  =  .01  with  to  =  5 , 
a  -  3 . 

(v)  A  Model  of  an  Overdense  Plasma. 

An  overdense  plasma  column  can  be  modelled  by  an  index  of 

refraction  which  becomes  negative,  such  as  the  example  used, 

2 

n  =  (r  -4)/5  for  r  <  3,  n  =  1  for  r  >  3.  Geometrical  optics 

2  2 

predicts  [9]  a  caustic  on  the  ellipse  ^  ^  =  1. 

Figure  4  shows  the  calculated  field  in  polar  coordinates 
with  w  =  5,  a  =  2,  At  =  0.05,  Ar  =  0.1,  and  A0  =  it/40.  Recalling 
that  the  larger  (hence  darker)  numbers  correspond  _o  constructive 
interference  note  that  the  ellipse  (heavy  curve)  lies  very  close 
to  the  edge  of  the  constructive  interference  for  tt/2  <  0  <  it  .  The 
rest  of  the  shadow  is  cast  by  the  caustic  along  the  line  y  =  r  sin0  =  3 
(remaining  portion  of  the  heavy  curve  for  r  >  3,  0  <  tt/2).  Figure 
5  shows  the  x-y  plot  of  the  interference  pattern  where  the  lighter 
regions  now  correspond  to  constructive  interference.  The  amplitude 
ranges  for  the  various  shades  are  again  (0.2),  ( 2 , 6 ) , ( 6 , 8) , ( 8 , 10 )  , 

( 10 , 14 ) ,  ( 14 , 17 )  ,  and  (17,°°).  Unfortunately  the  interpolation  re¬ 
quired  for  the  rectangular  output  smears  out  the  caustic. 


6.  Conclusions. 


There  are  four  obvious  dimensionless  parameters  in  the 
computed  problem,  wAr,  carA0,  wa  and  wr^  where  a  is  a 
relevant  length  and  r^  is  the  value  of  r  where  the  radiation 
condition  is  imposed.  The  sq”u  e  of  the  first  two  enter  the 
errors  produced  by  the  second  order  difference  scheme  and  in 
most  of  our  calculations  was  approximately  .2.  The  third 
measures  th~:  relevance  of  geometrical  optics  and  ranged  up 
to  15.  The  last  one,  usually  about  .  025  ,  arises  from  an 
expansion  at  00  for  fixed  u).  The  error  terms  are  of  order 
(wr^)  but  the  coefficients  are  singular  as  to 

In  many  of  the  problems  there  were  discontinuities, 
infinitely  thin  objects  or  discontinuities  in  the  index  of 
refraction.  These  induce  discontinuities  in  the  second 
derivatives  of  the  time  dependent  solution  and  in  the  steady 
state.  In  spite  of  this  source  of  inaccuracy  there  was 
remarkable  agreement  w^th  known  results  and  there  is  probably 
some  cancellation  of  error. 

The  method  could  be  applied  to  higher  dimensional 
phenomena  where  the  reflecting  object  as  well  as  variable 
index  of  refraction  are  not  amenable  to  solution  by  integral 
equations.  Since  the  problem  is  solved  by  a  very  straight 
forward  difference  scheme  it  should  be  easily  effected  by 
vector  computation.  Other  problems  where  variations  of  the 
limiting  amplitude  principle  could  be  used  involve  Maxwell's 
equations,  and  wave  propagation  in  water. 
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Finally,  the  method  has  been  modified  slightly  and  applied 
to  the  interesting  and  important  case  of  a  nonlinear  medium  [12]. 

In  that  report 

n  =  nQ  +  y (1  -  nQ)  |u|  2 

where  n^  is  the  index  of  refraction  used  in  our  overdense  plasma 
model  and  y  is  a  constant  which  controls  the  strength  of  the  non¬ 
linearity.  The  results  were  excellent;  the  effects  of  refraction 
and  self-focusing  could  be  discerned  and  controlled  by  varying 
y.  In  pai ticular,  self-focusing  amplifies  the  fields  near  the  origi 
as  was  first  observed  by  F.  Tappert  [unpublished  results]. 
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Figure  Captions 

Figure  1:  The  interference  pattern  produced  by  the  total  field's 

amplitude,  |U|,  for  a  metal  cylinder.  The  light  regions  correspond 
to  constructive  interference,  the  dark  to  destructive.  The  dark 
semicircle  is  the  cylinder.  The  relevant  parameters  are  toa  =  10, 

Ar  =  0.1,  AO  =  tt/40,  and  At  =  0.05. 

Figure  2 :  The  total  field's  amplitude,  |U|,  in  polar  coordinates 

for  a  metal  strip  with  the  incident  wave  normalized  to  10.  The  in¬ 
cident  wave  strikes  the  strip  head  on  (a  =  0).  The  set  of  points 
{ ( r ,  ©  )  |  0  =  it  /  2  ,  3tt/2,  0  <  r  <  2}  represents  the  strip  n  polar 
coordinates.  Since  the  numbers  printed  are  the  actual  field 
values,  the  dark  regions  now  correspond  to  constructive  inter¬ 
ference,  the  light  to  destructive.  The  relevant  parameters  are 
the  same  as  in  Figure  1  except  A0  =  tt/20. 

Figure  3:  The  interference  pattern  produced  by  the  total  field's 

amplitude,  |u|,  for  a  dielectric  lens  with  n  =  0.4  The  relevant 

parameters  are  ua  =  15,  Ar  =  0.1,  A9  =  tt/20,  and  At  =  0.05. 

Figure  4:  The  total  field's  amplitude,  |U|,  in  polar  coordinates 

for  an  overdense  plasma  column.  The  incident  wave  was  normalized 

2  2 

to  10.  The  heavy  curve  represents  the  caustic  ^  ^  =  1 ,  when 

r  <  3  and  the  geometric  shadow  y  =  3  when  r  >  3.  The  relevant 
parameters  are  ua  -  15,  Ar  =  0.1,  A0  =  tt/40,  and  At  =  0.05. 

Figure  5 :  Same  physical  problem  and  relevant  parameters  as 

Figure  4.  However',  this  is  now  the  interference  pattern  in  x-y 
coordinates  with  the  light  regions  corresponding  to  constructive 
interference,  the  dark  to  destructive. 
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Legend  for  Table  1:  Columns  A-F  present  the  cross  section 

£TUd)l/2  S  ( Q  )  ,  for  a  metal  cylinder  with  ug  ~  S  ( 0  )e^liir  ,  At  =  0.05, 
and  AG  =  tt/40.  /r 

Column  A  contains  the  results  given  in  [2],  B  presents  a 
tabulated  version  of  the  data  given  graphically  in  [1],  and 
columns  C-F  contain  the  results  of  our  ..->w  calculations. 

Legend  for  Table  2:  Columns  A-D  present  the  cross  section 

(Tro^l/2  2(0)^  for  a  metai  strip  with  ug  ~  S(6)e_^wr  ,  At  =  0.05, 

Ar  =  0.1,  and  AG  =  tt/20.  ^ 

Columns  A  and  D  present  tabulated  versions  of  data  given 
graphically  in  [1].  Columns  B  and  C  contain  the  results  of 
our  new  calculations.  The  parameter  a  is  the  angle  between 
the  incident  plane  wave  and  the  normal  to  the  strip.  L  is  the 
length  of  the  strip. 

Legend  for  Table  3:  Columns  A-C  present  the  cross  section 
S(G),  for  various  Helmholtz  resonators  with 

us  ~  S(0)ei“r  ,  At  =  0.05,  AG  =  tt/20,  and  Ar  =  0.1.  The 

parameter  a  is  the  angle  between  the  incident  plane  wave  and 
the  positive  x  axis  while  6  is  the  aperture  angle  of  the 
resonator.  The  parameter  wa  =  10  where  a  =  radius  of  the 


resonator . 
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